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Resumen 


Se aborda la descripción del movimiento del agua en un sistema de drenaje 
agrícola subterráneo con un modelo basado en el acoplamiento de las ecuaciones 
de Richards y de Boussinesq, con capacidad de almacenamiento variable sujeta 
en la frontera de los drenes a una condición de radiación fractal. El modelo se 
aplica para describir el funcionamiento hidráulico de un sistema de drenaje 
instalado en campo, el cual permite estudiar separaciones entre drenes de 40 y 
25 m. La mayor parte de la caracterización hidráulica del sistema es realizada de 
manera independiente de los eventos de drenaje, empleando una metodología 
basada en la curva granulométrica y porosidad total del suelo. Las evoluciones 
experimentales de la lámina drenada y del abatimiento del manto freático asociadas 
con la separación entre drenes de 40 m permiten calibrar dos parámetros del 
modelo propuesto, la escala de presión que interviene en la curva de retención de 
humedad del suelo y el parámetro de escala de la radiación fractal. La evidencia 
experimental correspondiente a la separación entre drenes de 25 m permite 
mostrar que el parámetro de escala de la radiación fractal es independiente de la 
separación entre drenes y que puede suponerse proporcional a la conductividad 
hidráulica a saturación del suelo. Se muestra que las hipótesis clásicas del drenaje 
agrícola que representan de manera simplificada las transferencias desarrolladas 
en la zona parcialmente saturada del suelo, considerando el término de recarga 
de la ecuación de Boussinesq nulo o igual a la evapotranspiración, trasladan la 
incertidumbre de lo que pasa en esta zona a los parámetros estimados a partir de 
las pruebas de drenaje. 


Palabras clave: zona vadosa, zona saturada, acuífero libre somero, recarga, 
descarga. 


Introducción 


El análisis detallado del flujo del agua en 
un sistema de drenaje agrícola subterráneo 
requiere considerar los procesos de transferencia 
desarrollados en las zonas parcialmente satu- 
rada y saturada del suelo. De estas zonas, es 
en la parcialmente saturada, también llamada 
zona vadosa, donde se desarrolla lo esencial de 
los procesos de infiltración, extracción de agua 


por las plantas, evaporación, almacenamiento y 
recarga de acuíferos. 

La mayor parte de los modelos hidrológicos 
que describen el flujo del agua a través del perfil 
del suelo sobre-simplifican la descripción de las 
transferencias desarrolladas en la zona vadosa 
(Kumar et al., 2008); estos modelos consideran 
relaciones de naturaleza esencialmente empí- 
rica para determinar el volumen de agua que 
se infiltra en el suelo, por ejemplo, el método 
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de capacidad de infiltración media, criterio 
del coeficiente de escurrimiento, criterio 
del Servicio de Conservación de Suelos de 
Estados Unidos y método de los números 
de escurrimiento (Chow et al., 1988). La 
naturaleza empírica de estos métodos origina 
una fuerte incertidumbre sobre la generalidad 
y el comportamiento de sus parámetros en 
condiciones distintas a las de su derivación. 

En un intento por mejorar la representación 
de la zona no saturada, otro grupo de modelos, 
como por ejemplo SWAT (Arnold et al., 1998), 
HEC-HMS (U.S. Army Corps of Engineers, 
2001) y SWMM 5.0 (Rossman, 2009), han 
incorporado relaciones con base teórica como 
las ecuaciones de Green y Ampt (1911), y 
Horton (1940), que describen la evolución en 
el tiempo de la lámina infiltrada; sin embargo, 
ambas relaciones consideran hipótesis que 
limitan su capacidad de representar los 
procesos desarrollados en un medio real. 

La descripción detallada del flujo del agua 
en las zonas parcialmente saturada y saturada 
del suelo puede hacerse con la ecuación de 
Richards (1931), que resulta de la combinación 
del principio general de conservación de 
masa y de la ley de Darcy generalizada. Sin 
embargo, la aplicación de esta ecuación a casos 
tridimensionales e incluso bidimensionales 
exhibe fuertes problemas, uno de éstos es 
la cantidad de información experimental 
requerida para determinar los parámetros de 
las características hidrodinámicas del suelo 
(curva de retención de humedad y curva de 
conductividad hidráulica). Otro problema 
es la solución numérica de la ecuación de 
Richards, ya que su convergencia y estabilidad 
requiere una discretización espacial fina del 
dominio de solución, un control riguroso del 
paso de tiempo y procedimientos iterativos 
que demandan gran cantidad de cálculos para 
linealizar y resolver el sistema de ecuaciones; 
esto se traduce en tiempos de cómputo que 
se incrementan exponencialmente conforme 
aumenta el tamaño del dominio de solución. 

Una manera de solventar los problemas de 
aplicabilidad de la ecuación de Richards en tres 


dimensiones consiste en utilizar la ecuación 
de Boussinesq de los acuíferos libres, que es 
resultado de integrar en la vertical la ecuación 
de Richards, asumiendo una distribución 
hidrostática de presiones (e.g. Bear, 1972). En 
esta ecuación, las transferencias desarrolladas 
en la zona parcialmente saturada del suelo son 
consideradas de manera simplificada a través 
de un coeficiente de almacenamiento. 

En varios modelos de drenaje agrícola, el 
coeficiente de almacenamiento de la ecuación 
de Boussinesq se ha supuesto constante 
(Dumm, 1954; McDonald y Harbaugh, 1988; 
Upadhyaya y Chauhan, 2000); sin embargo, esta 
hipótesis limita sensiblemente la descripción 
de las transferencias desarrolladas en la zona 
vadosa. Otros estudios han introducido un 
coeficiente variable que depende de la carga 
hidráulica, pero consideran relaciones de 
naturaleza empírica (e.g. Pikul et al., 1974; 
Pandey et al., 1992; Gupta et al., 1994; Samani et 
al., 2007). Recientemente, Hilberts et al. (2005) 
y Fuentes et al. (2009) han mostrado que el 
coeficiente de almacenamiento es función de la 
carga hidráulica y han establecido su relación 
con la curva de retención de humedad. 

Cuando la recarga o descarga vertical al 
acuífero es nula, considerar el coeficiente de 
almacenamiento variable permite obtener 
representaciones aproximadas de la dinámica 
del agua en un acuífero libre somero (Fuentes 
et al., 2009); sin embargo, cuando la recarga o 
descarga es variable, las transferencias que 
se desarrollan en la zona vadosa son más 
complejas y se requiere considerar la variación 
de la conductividad hidráulica que interviene 
en la ley de Darcy generalizada; esto puede ser 
realizado acoplando la ecuación de Richards 
para la zona parcialmente saturada con la 
ecuación de Boussinesq para la zona saturada. 
de Richards 
proporciona la intensidad de la recarga o 


La solución de la ecuación 


descarga del acuífero, que es una variable 
de la ecuación de Boussinesq, mientras que 
la solución de la ecuación de Boussinesq 
proporciona la evolución de la superficie 
libre, que define las alturas de las columnas 
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de suelo parcialmente saturado donde se 
aplica la ecuación de Richards. La hipótesis en 
este acoplamiento es que el movimiento en la 
zona vadosa es esencialmente vertical, lo cual 
permite considerar la forma unidimensional 
de la ecuación de Richards, y evitar cálculos y 
tiempos de cómputo excesivos. 

Pikul et al. (1974) desarrollaron un modelo 
para drenaje agrícola que acopla las ecuaciones 
de Richards y Boussinesq en su forma 
unidimensional, sin embargo consideran una 
relación para el coeficiente de almacenamiento 
empírica y difícil de determinar (Vachaud 
y Vauclin, 1975), y además imponen en los 
drenes condiciones de frontera simplificadas 
(tipo Dirichlet o carga hidráulica especificada). 
García et al. (1995) desarrollaron un modelo 
para describir la dinámica del agua y el 
transporte de solutos en medios de saturación 
variable, acoplando las ecuaciones de Richards 
unidimensional, Boussinesq y advección-dis- 
persión; recientemente, Kumar et al. (2008) han 
acoplado el modelo HIDRUS 1D (Simúnek et al., 
2005), que resuelve la ecuación unidimensional 
de Richards, con el modelo MODFLOW-2000 
(McDonald y Harbaugh, 1988), que resuelve la 
ecuación general de los acuíferos y en particu- 
lar la ecuación de Boussinesq para acuíferos 
libres someros. Sin embargo, en ambos aco- 
plamientos se tiene como limitante que el 
coeficiente de almacenamiento es considerado 
constante y que la condición de frontera 
usada en los drenes agrícolas es del tipo 
Dirichlet (carga hidráulica especificada) o de 
radiación lineal (flujo de drenaje proporcional 
a la diferencia de presiones en la vecindad del 
dren). 

En relación con la descripción del flujo de 
drenaje, se han realizado diversos estudios 
sobre el efecto del tamaño, forma y distribución 
de las perforaciones de la pared del dren en 
la formación de una carga hidráulica sobre 
éste y en la evacuación de agua del suelo 
(Kirkham, 1950; Engelund, 1953; Ernst, 1962; 
Youngs, 1974; Skaggs, 1978). Algunos modelos 
de la literatura relacionan las propiedades 
del sistema suelo-dren con la resistencia al 


flujo del agua en su interfaz, y establecen 
relaciones lineales o polinomios de segundo 
grado entre el gasto drenado y la carga 
hidráulica (Kohler et al., 2001). Zavala et al. 
(2007) analizan el comportamiento del flujo de 
drenaje, considerando el suelo y la pared del 
dren como objetos fractales, y establecen que la 
transferencia de agua del suelo al interior del 
dren debe ser representada con una condición 
de frontera de radiación no lineal o fractal. En su 
análisis consideran la ecuación de Boussinesq 
con coeficiente de almacenamiento constante, 
recarga nula, y calibran los parámetros de 
la radiación fractal a partir de un evento de 
drenaje de laboratorio. 

El objetivo de este trabajo es desarrollar un 
modelo para el drenaje agrícola subterráneo 
que considere las transferencias desarrolladas 
en la zona vadosa del suelo bajo recarga 
o descarga vertical variable en el tiempo, 
basado en el acoplamiento de la ecuación de 
Richards unidimensional y la ecuación de 
Boussinesq unidimensional con coeficiente 
de almacenamiento variable sujeta en los 
drenes a una condición de radiación fractal. Se 
realizará la estimación de los parámetros de la 
radiación fractal considerando datos obtenidos 
en una prueba de drenaje realizada en campo 
y se evaluará su capacidad de descripción 
considerando datos medidos en un segundo 
evento de drenaje. 


Ecuaciones de base 


La figura 1 muestra un esquema del flujo no 
saturado-saturado que se presenta en un 
sistema de drenaje subsuperficial; se considera 
un acuífero libre somero con drenes paralelos 
equidistantes y estrato impermeable horizon- 
tal. En esta figura se representa una columna 
de suelo parcialmente saturado, cuya altura 
varía en el tiempo y espacio, en función de las 
condiciones de flujo en la superficie del suelo, 
del flujo de agua que se evacua a través de los 
drenes, de las propiedades físicas e hidráulicas 
del medio poroso y de las características 
geométricas del sistema. La base de cada 
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columna de suelo corresponde a la posición de 
la superficie libre de agua o nivel freático, y su 
altura se extiende hasta la superficie del suelo. 


Zona parcialmente saturada (zona vadosa) 


La descripción del movimiento del agua en 
una columna de suelo parcialmente saturado, 
homogéneo e isotrópico, puede realizarse con 
la ecuación de Richards (1931), que resulta de 
la combinación de la ecuación de continuidad 
y la ley de Darcy generalizada: 


Kw (5 +1) a) 


donde y = vy(z, t) es el potencial de presión 
del agua en el suelo, expresado como una 
altura equivalente de columna de agua [L]; el 
potencial gravitacional ha sido asimilado a la 
coordenada vertical z, orientada positivamente 
en dirección ascendente [L]; t es el tiempo [T]; 
K(w) es la conductividad hidráulica en fun- 


C(w) =d0(w) / dy esla capacidad específica [L”]; 
9 = O(y) es el contenido volumétrico de agua o 
contenido de humedad [L*L”] y depende de la 
presión del agua en el suelo, esta dependencia 
es conocida como la característica de humedad 
o curva de retención de humedad. 


Condiciones límite 


La resolución de la ecuación de Richards 
requiere la definición de las condiciones límite 
que representen el proceso que se desarrolla 
en la columna de suelo parcialmente saturado. 
Como condición inicial debe especificarse la 
distribución de presiones en la zona vadosa del 
suelo: 


y(2,0) =y(z) (2) 


La base de la columna de suelo correspon- 
de a la posición de la superficie libre de agua 
(potencial de presión nulo), por lo que debe 
imponerse una condición de frontera tipo 
Dirichlet homogénea: 


ción de la presión del agua en el suelo [LT]; y=0 z=H y 1t>0 (3) 
z 
Columna de suelo 
Superficie 
Hs 
Zona parcialmente saturada 
Nivel freático 

EN A AAA Dren 


Zona saturada 


Estrato impermeable 


Figura 1. Esquema del flujo parcialmente saturado-saturado en un sistema de drenaje. 
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Hes la elevación de la superficie libre o nivel 
freático, medida desde el estrato impermeable. 
En la frontera superior de la columna de 
suelo puede usarse una condición tipo Dirichlet 
si se desea representar el efecto de saturación 
asociado con el riego por gravedad o al 
generado por un escurrimiento superficial: 
y=0 z=H, y t>0 (4.1) 
O bien puede usarse una condición de 
frontera tipo Neumann de flujo prescrito: 


a. 40) z=H, y t>0(42) 


donde H, es la elevación de la superficie del 
suelo medida desde un nivel de referencia, en 
este caso el estrato impermeable; y f (t) es un 
valor de flujo conocido, como la intensidad 
de la evaporación o de la precipitación. La 
condición (4.2) también puede usarse para 
representar la evapotranspiración de un 
cultivo, sobre todo en la etapa inicial de su 
ciclo fenológico, cuando el desarrollo de las 
raíces no es profundo y la mayor parte del 
proceso evapotranspirativo se desarrolla en 
zonas del suelo cercanas a la superficie. 


Propiedades hidráulicas del suelo 


La resolución de la ecuación de Richards 
(1) sujeta a las condiciones (2-4) requiere el 
uso de modelos que relacionen el contenido 
volumétrico de agua 8 con el potencial de 
presión del agua en el suelo y, así como la 
conductividad hidráulica K con el contenido 
volumétrico de humedad 6. Fuentes et al. 
(1992) han analizado los modelos clásicos 
usados en la literatura para representar las 
propiedades del suelo y han demostrado que 
la mayor parte no satisface las propiedades 
integrales de la infiltración; también han 
determinado que una combinación de modelos 
que satisface estas propiedades y que es útil en 
estudios experimentales es el modelo para la 


característica de humedad de van Genuchten 
(1980), sujeto a la restricción de Burdine (1953): 


o(y)-0, _ y) 


9, = O, Va 


=m 


m=1-2/n (5) 


donde 6, es el contenido volumétrico de agua a 
saturación [L*L*]; 0 es el contenido volumétrico 
de agua residual [L*L”]; y, es un parámetro de 
escala de la presión [L]; m y n son parámetros 
de forma empíricos adimensionales. 

Y el modelo para la conductividad hidráu- 
lica de Brooks y Corey (1964): 


K0)=x| cl ] (6) 


9, -0, 


donde K._ es la conductividad hidráulica a 
saturación [LT*] y n es un parámetro forma 
adimensional. 


Zona saturada 


La dinámica del agua en un acuífero libre 
somero puede ser representada con la ecuación 
de Boussinesq del drenaje agrícola, que resulta 
del principio de conservación de masa, la ley 
de Darcy y la hipótesis de Dupuit-Forcheimer 
concerniente a una distribución hidrostática de 
presiones: 


o 


donde T = K(H — H) es la transmisibilidad del 
acuífero [L*T*]; (H — H,) representa el espesor 
del acuífero [L], siendo H y H, las elevaciones 
de la superficie libre o carga hidráulica y del 
estrato impermeable [L], medidas desde el 
mismo nivel de referencia; cuando el estrato 
impermeable es aproximadamente horizontal, 
se puede asumir H, = 0; R es el volumen de 
recarga O descarga vertical por unidad de área 
de acuífero por unidad de tiempo [L*L?T"); 
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y u(H) es el coeficiente de almacenamiento 
[E*L?], que en un acuífero libre es función de la 
carga hidráulica (Hilberts et al., 2005; Fuentes 
et al., 2009). 


Coeficiente de almacenamiento 


De acuerdo con Fuentes et al. (2009), el 
coeficiente de almacenamiento en un acuífero 
libre somero está relacionado con la curva de 
retención de humedad de la siguiente manera: 


0.-0(H-H,) ¿H_<H,, 


8 
0,- 0 (HH) + (0,-0/.); Hg <H, / 


ue 


donde HA, = H + |w,| y y, es el potencial de 
presión asociado con el concepto de capaci- 
dad de campo, en suelos de textura arenosa 
|y,] =1.0 m y suelos arcillosos |y,| = 3.3 m. 

La introducción de la curva de retención 
(ecuación (5)) en la definición del coeficiente de 
almacenamiento (ecuación (8)) permite obtener 
la siguiente relación: 


0-0) (005) val 7") 5,55, 
E 
+(0, 0); H¿<H, 
Condiciones límite 


De acuerdo con Zavala et al. (2007), la descrip- 
ción de la dinámica del agua en sistemas 
de drenaje subterráneos con la ecuación (7) 


5 requiere el uso de las siguientes condiciones 
» límite: 

2 

a H(:0)=H,() (10) 
E A, z=0 y x=E £>0 (1) 
dE dx 

S, donde H,(x) es la posición inicial del manto 


freático, que puede ser función de la coordenada 
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horizontal x. El signo positivo en la condición 
de frontera (11) se toma para el dren ubicado en 
x = 0, mientras que el negativo para el que está 
ubicado en x = L, donde L es la separación entre 
drenes. En la ecuación (11) se introduce el flujo 
de drenaje (q,) proporcionado por la radiación 
fractal: 


S ym ES, 
H-D m*sd 
=Q0| — 12 
Ya 0 P, ] ( ) 


q, es el valor máximo del flujo de drenaje; 
P,es la profundidad de los drenes; D, es la 
profundidad del estrato impermeable medida 
desde la posición de los drenes; s, y s, son la 
dimensión cociente del suelo y de la pared del 
dren (razón entre la dimensión fractal del objeto 
D, y la dimensión del espacio de Euclides). La 
dimensión cociente de un objeto fractal (s) está 
relacionada con su porosidad volumétrica (6) 
de acuerdo con: 


(1-4 +09 =1, s=D,/3 (13.1) 


y considerando que la porosidad areal es 
9 = 4”, se obtiene: 


EN 
(= o) +Q2=1 (13.2) 


En las ecuaciones (13.1) y (13.2) se 
introduce s,, y s, en lugar de s para calcular las 
dimensiones del suelo y de la pared del dren, 
respectivamente. 

El gasto total evacuado por cada dren es 
proporcionado por la relación: 


Q =2/(D, +H)g, (14) 
donde l es la longitud del dren. 
Solución numérica 
Los sistemas de ecuaciones (1-6) y (7-14) se 
resuelven numéricamente de la siguiente 
manera: la discretización espacial se realiza 


con el método del elemento finito tipo Galerkin 
(Zienkiewicz et al., 2005); la discretización 
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temporal se hace aplicando un método de 
diferencias finitas (Herbert y Anderson, 1995); 
los sistemas resultantes de la discretización 
se linealizan empleando el procedimiento 
iterativo de Picard, mientras que los sistemas 
resuelven 
usando el método de gradiente conjugado 


de ecuaciones algebraicas se 
precondicionado propuesto por Noor y Peters 
(1987). 

El esquema numérico para la ecuación de 
Richards es: 


Act p p+l 
At + —B HOAt Wions 
(15) 


AP 
== os == ( E, tv oi oa 


donde t es el tiempo; At es el paso de tiempo; 
es un coeficiente de ponderación en el tiempo 
(0 < O < 1); p indica el número de iteración en 
el intervalo de tiempo; A es la matriz de masa; 
B es la matriz de rigidez; G es el vector de la 
fuerza de gravedad, F es el vector de flujos en la 
frontera (z = H.); y es el vector que contiene los 
valores del potencial de presión en los nodos 
del elemento finito. El uso de funciones de 
interpolación lineales y funciones del sistema 
de masa concentrado (Zienkiewicz et al., 2005) 
permite obtener los coeficientes siguientes: 


Aj=8,Y.C) > (16.1) 
B, = zon” : (16.2) 
as 2 Mm (16.3) 
E=tjls ¿=9 (2=4H) (16.4) 


n es el número de nudos considerados en la 
discretización de la columna de suelo; e = n — 
1 es el número total de elementos finitos; i y j 
son índices para referenciar el número de nudo 
de cada elemento (1, 2,..., 1), el primer índice 


está asociado con la solución aproximada 
de elemento finito, mientras que el segundo 
índice está relacionado con las funciones de 
peso que se utilizan para minimizar el error 
de la discretización espacial (Zienkiewicz 
et al., 2005); Ó, es la delta de Kronecker; C es 
la capacidad específica en los nudos del 
elemento; Az es el tamaño del elemento finito; 
Kesla conductividad hidráulica en el elemento, 
tomada como el promedio aritmético de las 
conductividades en los nudos del elemento); f (t) 
es la evolución del flujo de Darcy en la superfi- 
cie del suelo; el signo positivo se toma cuando 
se representa la intensidad de la precipitación 
y el signo negativo cuando se representa la 
evaporación o evapotranspiración en función 
del tiempo. 

El esquema numérico para la ecuación de 
Boussinesq es el siguiente: 


AP 1 
ES Ea Bo 
(17) 
A p 
5 e -(1—00)B sy et E+OAt 


donde H es el vector que contiene los valores de 
la carga hidráulica. En este caso, los coeficientes 
de la matriz de masa M, matriz de rigidez B, 
vector de flujos en la frontera F y vector de 
recarga o descarga R_ son: 


A 
A¡=5,Y hy 3 (18.1) 
B.= Nay E (18.2) 
U) - Ax 
G:==q4(H)l: ¿=1, i=N (18.3) 
Ax 
R.= NR 18.4 
RS (18.4) 


e 


N es el número de nudos considerados 
en la discretización del dominio de solución 
definido a partir de la separación entre drenes 
L; 1 y ¡ son índices para referenciar el número 
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de nudo del elemento finito (1, 2,..., N), el 
primer índice está asociado con la solución 
de elemento finito, mientras que el segundo 
índice está relacionado con las funciones 
de peso que se utilizan para minimizar el 
error de la discretización espacial; e = N — 1 
es el número de elementos finitos: Ax es el 
tamaño del elemento; u es el coeficiente de 
almacenamiento en los nudos del elemento; T es 
la transmisibilidad en el elemento, tomada como 
el promedio aritmético de las transmisibilida- 
des en los nudos del elemento finito; q, es el 
flujo de drenaje calculado con la radiación 
fractal (ecuación (12)); y R es la intensidad de 
la recarga o descarga del acuífero, que es igual 
a la magnitud del flujo de Darcy sobre la base 
de la columna de suelo parcialmente saturado 


q = H,t) =-K(wlloy/oz) + 1]. 
Acoplamiento numérico 


a) En el intervalo de tiempo At se resuelve 
primeramente el sistema asociado con la 
ecuación de Boussinesq, considerando un 
valor de recarga o descarga (R): cero si es el 
tiempo inicial (t = Af) y la primera iteración 
(p = 1); el valor calculado en el tiempo 
anterior sit>At y p= 1; o el valor obtenido 
en la iteración previa si p > 1. Con este 
procedimiento se obtiene la posición de la 
superficie libre en el intervalo de tiempo. 

b) A partir de la posición del manto freático 
H(x,t) se definen las alturas de las columnas 
de suelo parcialmente saturado y se 
resuelve el sistema asociado con la ecuación 
de Richards en cada una de estas columnas, 
obteniéndose la distribución de presiones y 
la magnitud del flujo de Darcy en la base de 
las columnas. 


Los pasos a y b se repiten en el intervalo 
de tiempo At hasta satisfacer el criterio de 
convergencia seleccionado (tolerancia para 
el máximo cambio tanto de la posición de la 
superficie libre como del flujo de Darcy en la 
base de las columnas de suelo entre iteraciones 


consecutivas); en este proceso iterativo se 
ajusta continuamente la altura de las columnas 
de suelo parcialmente saturado en función de 
la posición de la superficie libre. 


Aplicación 


El modelo Richards-Boussinesq se aplica para 
caracterizar y describir el funcionamiento del 
sistema de drenaje subterráneo de una parcela 
agrícola localizada en el municipio de Sinaloa 
de Leyva, Sinaloa, México, perteneciente al 
distrito de riego 075 Río Fuerte. La parcela 
tiene una superficie de 8.4 ha y cuenta con 
seis líneas paralelas de drenes subterráneos de 
polietileno corrugado de 10 cm de diámetro, 
instaladas a una profundidad media P, = 1.20 
m, que descargan libremente a un dren colector 
superficial; tres de las líneas presentan una 
separación entre drenes L, = 40 m y drenan el 
68% de la superficie total, mientras que las tres 
restantes tienen una separación entre drenes 
L, =25 m y cubren el 32% de la parcela. La pared 
de los drenes tiene entre cada corrugación seis 
perforaciones distribuidas a la largo de su perí- 
metro (cada 60”); el área perforada por metro 
de dren es de 32 cm”, por lo que la porosidad 
areal de la pared es d, = 0.0102 cm? /cm?. 

Para obtener información experimental 
que permita aplicar el modelo numérico, se 
monitoreó el funcionamiento de las líneas de 
drenaje centrales asociadas con cada separa- 
ción entre drenes, lo cual elimina las pertur- 
bacionesoriginadas porelcambio deseparación. 
Utilizando el plano de la parcela experimental, 
se delimitó el área de influencia de cada una 
de las líneas de drenaje seleccionadas, se 
dividió cada zona en tres secciones y en cada 
celda se seleccionó un punto de muestreo de 
suelo, considerando tres profundidades: 0-35 
cm, 35-70 cm y 70-130 cm. En el perfil de suelo 
se detectó un estrato de suelo prácticamente 
impermeable a 4.20 m de profundidad, por lo 
que D,=3.0 m. 

Las muestras de suelo, nueve por cada 
línea de drenaje, se analizaron en laboratorio 
para determinar la porosidad, granulometría 
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y textura; los resultados obtenidos muestran 
que se tiene poca variabilidad espacial de 
estas propiedades físicas en cada una de las 
zonas analizadas (cuadro 1). En consecuencia, 
se justifica considerar en la modelación la 
hipótesis de suelo homogéneo para cada zona 
de suelo. 

Previo a la evaluación hidráulica del 
sistema de drenaje se realizaron las siguientes 
actividades: rastreo, trazo de regaderas y 
colectores siembra en seco, 
trazo de surcos a 0.80 m de separación y 


superficiales, 


construcción de pozos de observación del 
manto freático a la mitad de separación 
entre drenes. Posteriormente se taparon las 
descargas de las seis líneas de drenaje y se 
aplicó el riego superficial de germinación hasta 
saturar completamente el suelo. Concluida la 
aplicación del riego se mantuvieron cerradas 
las descargas de las líneas de drenaje durante 
doce horas para permitir una mejor distribución 
del agua y de la presión en el suelo. Finalmente 
se destaparon las seis líneas de drenaje y se 
midió durante 15 días el volumen de agua 
evacuado a través de los drenes centrales 
asociados con la separaciones de 40 y 25 m; de 
manera simultánea se midió la evolución del 
manto freático a la mitad de separación entre 
drenes. 

La evolución en el tiempo de la 
evapotranspiración de referencia (Et) en 
el periodo del experimento se tomó de los 
registros de la estación climatológica automá- 
tica Batequis, localizada a 6 km de la parcela 


experimental. Esta información sirvió de base 
para estimar la evapotranspiración del trigo 
sembrado en la parcela experimental (Et ); el 
valor del coeficiente de cultivo considerado en 
la estimación K, = 0.40 corresponde a la etapa 
inicial del trigo. En la figura 2 se muestra la 
variación en el tiempo de la evapotranspiración 
de referencia y la evapotranspiración del trigo. 


Caso 1. Separación entre drenes de 40 m 


Se aborda la descripción de los procesos de 
transferencia desarrollados en la zona de suelo 
dominada por la línea de drenaje central de la 
separación entre drenes L, = 40 m. La longitud 
total de este dren es / = 446.20 m y su pendiente 
longitudinal es S, = 0.12%. 

La modelación del evento de drenaje 
requiere la determinación de los parámetros 
que intervienen en la curva de retención de 
humedad (ecuación (5)), curva de conduc- 
tividad hidráulica (ecuación (6)), coeficiente 
de almacenamiento variable (9) y condición de 
radiación fractal (12). Es esencial minimizar el 
número de parámetros a determinar a partir 
del evento transitorio de drenaje monitoreado, 
ya que entre mayor es la cantidad de paráme- 
tros que se estimen de una sola prueba, se gana 
en cuanto a la descripción de la experiencia 
particular, pero se pierde en cuanto al sentido 
de generalidad y capacidad de explicación 
del porqué y del cómo de los procesos. En 
consecuencia, la caracterización hidrodinámica 
se realiza aplicando la metodología reco- 


Cuadro 1. Valores medios de las propiedades físicas del suelo de la parcela experimental. = 

A 

”"u 

uv 

P i Arcill Li A 5 

dates io rcilla imo rena SEN 2 
cm'/cm % % 7 EE 

a 

Nueve muestras de suelo representativas de la zona con separación L = 40 m y 

Valor medio 0.469 42.60 34.60 22.80 Arcillosa ES 

em 

Desviación estándar 0.014 4.25 5.17 5.40 - E 
Nueve muestras de suelo representativas de la zona con separación L = 25 m = 

9 

Valor medio 0.474 45.70 30.85 23.45 Arcillosa e 
50 

Desviación estándar 0.007 4.90 1.90 3.70 = = 
== 
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Evapotranspiración (mm /d) 


A 


Figura 2. Evapotranspiración de referencia Ef, (estación climatológica Batequis) y evapotranspiración del trigo Et. 


(parcela experimental) estimada con un coeficiente de cultivo K.=0.4. 


mendada por Fuentes (1992), basada en la curva 
granulométrica y porosidad total del suelo. 

El contenido volumétrico de agua a 
saturación 6, se asimila a la porosidad total del 
suelo $. De acuerdo con los datos del cuadro 
1 se tiene una porosidad promedio para 
esta zona de suelo de q = 0.469 cm*/cm?, en 
consecuencia, 0, = 0.469 cm?/cm?. De acuerdo 
con Haverkamp et al. (2005), para reducir los 
parámetros de ajuste, el contenido volumétrico 
de agua residual puede asumirse igual a cero 
(9, = 0 cm*/cm?), lo cual implica aceptar que 
cuando y > —« necesariamente 9 > 0. 

El parámetro m de las relaciones (5) y (9), 
así como el parámetro 1 de la relación (6), se 
estiman a partir de la curva granulométrica 
del suelo y de su porosidad total. De acuerdo 
con el modelo funcional de la curva de reten- 
ción (5), se ajusta la curva granulométrica del 
suelo con la función de distribución: 


-M 


DY 
F(D)= m7) (19) 


donde F(D) es la frecuencia acumulada basada 
en el peso de las partículas, cuyos diámetros 
son inferiores o iguales a D; D, es un tamaño 
característico de las partículas; M y N son 
dos parámetros de forma empíricos, que de 
acuerdo con la restricción de Burdine (1953) se 
relacionan como M=1-2/N, con0<M=<1 
y N >2. El mejor ajuste de la ecuación (19) con 
la curva granulométrica representativa de 
la zona en análisis (R? = 0.9975) se obtiene con 
D, = 280 um y M =0.084 (figura 3). 

Haciendo A = mn y m, = MN, el parámetro 
de forma m de la curva de retención de van 
Genuchten puede relacionarse con M a través 
de la fórmula m,/4 = 1 + (Qs, - 1)/2(1 - s,). 
Introduciendo el valor de la porosidad volu- 
métrica del suelo en la ecuación (13.1) se tiene 
s,, = 0.688, y con el valor de m,= MN = 0.183 
se deduce 4 = mn = 0.115, de donde m = 0.054. 
La potencia de la curva de conductividad 
hidráulica se calcula con la relación y = 2s,, 
(1 + 2/12), obteniéndose n = 25.4. 

La conductividad hidráulica a saturación 
(K.) se determinó en laboratorio mediante una 
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Figura 3. Ajuste de la curva granulométrica de una muestra representativa de suelo de la zona 
L, = 40 m con la ecuación (19): D, =280 un y M= 0.084 (R? = 0.9975). 


prueba de carga de agua constante aplicada 
en una columna de suelo saturado; el valor 
obtenido es K_= 0.060 m/d. La porosidad 
volumétrica de la pared del dren se determina 
introduciendo su valor de porosidad areal en 
la relación (13.2) y resolviendo esta función 
implícita, obteniéndose s,= 0.569. 

Los parámetros que 
caracterización del sistema de drenaje son 
y, (ecuación (5)) y q, (ecuación (12)), y se 
determinan a partir de la evolución de la 


completan la 


lámina drenada experimental y la evolución 
del abatimiento del manto freático medida 
a la mitad de separación entre drenes. El 
procedimiento consiste en aplicar el modelo 
numérico Richards-Boussinesq (ecuaciones 
(15) a (18)) y determinar la combinación de 
valores que minimiza la diferencia entre datos 
medidos y calculados. En la modelación, 
la  evapotranspiración se 
manera 


representa de 
simplificada, imponiendo en la 
frontera superior de las columnas de suelo 
parcialmente saturado (superficie del suelo), 


una condición de frontera tipo Neumann 


variable en el tiempo (ecuación (4.2)), cuyos 
valores corresponden a la evapotranspiración 
del trigo estimada previamente (Et en la 
figura 2). La condición inicial de la simulación 
corresponde a la observada en el experimen- 
to, H(x,0) =H.. 

La combinación de valores que minimiza 
la diferencia entre las evoluciones de la lámina 
drenada y la superficie libre, experimental y 
calculada, es y, = 0.70 m y q, = 27.5 K, siendo 
el error cuadrático medio RECM = 1.62 x 10*m 
y RECM = 4.02 x 10* m, respectivamente. En 
las figuras 4a y 4b se muestran los resultados 
obtenidos. La correspondencia entre datos 
medidos y observados es satisfactoria. 

Para mostrar la importancia de la zona 
vadosa en el funcionamiento de un sistema 
de drenaje, se comparan los resultados del 
modelo contra los 
resultados del modelo basado sólo en la 
ecuación de Boussinesq con coeficiente de 
almacenamiento variable; se 


Richards-Boussinesq 


consideran 
los parámetros del suelo y de la radiación 
fractal obtenidos precedentemente, así como 
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Lámina de agua drenada (cm) 


O Experimental 


0 5 10 15 0 5 10 15 
Tiempo (día) Tiempo (día) 


Profundidad de la superficie libre (m) 


a) Evolución de la lámina drenada (RECM = 1.62 x 10% m) b) Evolución del manto freático a la mitad de separación 
entre drenes (RECM = 4.02 x 10% m) 


Figura 4. Calibración de los parámetros de escala de la curva de retención de humedad y de la radiación fractal: 
y,=-0.70 m y q,=27.5 K.. 


la condición de carga hidráulica observada Boussinesq, los valores que se obtienen en la 
al inicio del experimento. Se analizan dos calibración absorben la sobre-simplificación 
escenarios clásicos considerados en la lite- de los procesos de transferencia desarrollados 
ratura cuando se modela con la ecuación de en la zona parcialmente saturada del suelo; 
Boussinesq; el primero consiste en asumir que en este caso, los parámetros y, y q, hubiesen 
la descarga del acuífero somero por efecto tenido un valor diferente al determinado con 
de la evapotranspiración es nula (R = 0), el modelo Richards-Boussinesq. 

y el segundo consiste en aceptar que la eva- 

potranspiración del cultivo tiene un efecto Caso 2. Separación entre drenes de 25 m 
instantáneo sobre el manto freático, sin 

importar las condiciones de humedad, La segunda configuración del sistema de 
capacidad de retención y flujo en la zona drenaje experimental, L, = 25 m, es analizada 
vadosa, es decir R = Et (t). Los resultados de con el modelo Richards-Boussinesq. La línea 
las modelaciones numéricas se muestran en evaluada tiene una longitud de 353.50 m y 
las figuras 5a y 5b, donde se observa que el una pendiente de 0.12%; la profundidad 
modelo basado en la ecuación de Boussinesq del estrato impermeable medida desde 
con la hipótesis de descarga nula sobre-predice los drenes es D, = 3.0 m, y la variación en 
la evolución la lámina de agua drenada y el tiempo de la evapotranspiración usada 
subestima el abatimiento del manto freático; como condición de frontera superior de las 
mientras que con el supuesto de descarga igual columnas de suelo no saturadas es la que se 
a la evapotranspiración variable en el tiempo muestra en la figura 2. 

se subestima sensiblemente la evolución La caracterización hidrodinámica del sue- 
de la lámina drenada y se sobreestima el lo se realiza con el procedimiento descrito 
abatimiento del nivel freático. Estos resultados previamente. Los resultados obtenidos son 
permiten mostrar que cuando se realiza la 4 = 0.474 cm? / cmd, 9, = 0.474 cm* / cm, 9,=0.0 
estimación de parámetros de un sistema de cm?/cm?, m = 0.043, n = 32.25 y K.=0.082 m 1d. 
drenaje exclusivamente con la ecuación de Mientras que los parámetros de forma de la 
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b) Evolución de la superfice libre a la mitad de separacion 
entre drenes (L¡/2 =20 m) 


Figura 5. Comparación de los resultados obtenidos con el modelo Richards-Boussinesq y el modelo basado sólo 


en la ecuación de Boussinesq con descarga nula R = 0 y descarga igual a la evapotranspiración R =-Et. (t). 


radiación fractal son s, = 0.691, s, = 0.569. Se 
retiene con fines de validación el parámetro 
de escala de la condición de radiación fractal 
previamente estimado (q,=27.5 K). 

En esta modelación, el único parámetro 
de calibración es y, el cual se determina a 
partir de la evolución de la lámina drenada 
experimental. El valor que minimiza la 
diferencia entre la lámina experimental medida 
y la lámina calculada con el modelo Richards- 
Boussinesq es y, = 0.49 m (figura 6a), siendo 
la suma del error cuadrático medio RECM 
= 7.42 x 10* m. Los resultados presentados 
en la figura 6b muestran que la descripción 
del abatimiento del nivel freático obtenida 
con el modelo Richards-Boussinesq describe 
aceptablemente la evolución experimental 
(RECM = 5.09 x 10* m), sobre todo si se 
considera que la evolución medida no es 
elemento de calibración sino de validación; 
esta evidencia muestra que para el intervalo 
de separación entre drenes considerado en 
este trabajo (25 m < L< 40 m), el parámetro de 
escala de la radiación fractal q, no exhibe una 
dependencia aparente de la separación entre 


drenes y puede ser relacionado con el valor de 
K. 


Ej 


Conclusiones 


Los procesos de transferencia desarrollados 
en la zona parcialmente saturada del suelo o 
zona vadosa son de fundamental importancia 
para definir la dinámica del agua en el 
espesor saturado del suelo. En el caso de los 
sistemas de drenaje agrícola subterráneo, 
donde el abatimiento y/o control del manto 
freático somero es fundamental, es necesario 
representar estos procesos de manera explí- 
cita en los modelos matemáticos que describen 
su funcionamiento hidráulico. 

A fin de considerar los procesos desa- 
rrollados en la zona no saturada del suelo en 
el análisis del drenaje agrícola subterráneo, 
se ha desarrollado un modelo de simulación 
basado en el acoplamiento numérico de las 
ecuaciones de Richards y de Boussinesq 
con capacidad de almacenamiento variable 
sujeta en los drenes a una condición de 
radiación fractal. Se ha implementado un 
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Figura 6a. Calibración de la lámina drenada medida y la 
lámina descrita con el modelo Richards-Boussinesq 
con y, = 0.49 m (RECM = 7.42 x 10* m). 
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Figura 6b. Comparación de la evolución del manto freático 
a la mitad de separación entre drenes (L,/2 =12.5 m) 
observada y descrita con el modelo Richards-Boussinesq. 


algoritmo iterativo para la convergencia del 
acoplamiento Richards-Boussinesq, en el 
que se calcula y ajusta sistemáticamente la 
posición de la superficie libre, la altura de las 
columnas de suelo parcialmente saturado y 
la magnitud del flujo de Darcy en la base de 
estas columnas. 

El modelo desarrollado ha sido aplicado 
para caracterizar y describir el funcionamiento 


hidráulico del sistema de drenaje subterráneo 


de una parcela agrícola experimental del 
distrito de riego 075 Río Fuerte, Sinaloa, los 
drenes en esta parcela están instalados a una 
profundidad media de 1.20 m y presentan 
separaciones de 40 y 25 m, siendo el suelo de 
textura esencialmente arcillosa. Se hicieron 
dos pruebas de drenaje para determinar los 
parámetros del modelo y validar parcialmente 
sus resultados. 

Los parámetros de forma de las característi- 
cas hidrodinámicas del suelo y de la condición 
de radiación fractal usada en los drenes se han 
estimado empleando una metodología basada 
en la curva granulométrica y porosidad total 
del suelo. A partir de los datos de la prueba de 
drenaje asociada con la separación entre drenes 
L, = 40 m, se han determinado los valores de 
los parámetros de escala que intervienen en 
la curva de retención de humedad del suelo 
y en la condición de radiación fractal; el 
procedimiento consistió en minimizar el error 
entre las evoluciones medidas de la lámina 
drenada y la superficie libre, y las evoluciones 
calculadas con el modelo Richards-Bous- 
sinesq. La no presencia de oscilaciones e 
inestabilidades en las evoluciones de lámina 
drenada y de la superficie libre descritas con 
el modelo es un indicador de la consistencia 
del método numérico empleado para resolver 
el acoplamiento Richards-Boussinesq. 

La modelación de la prueba de drenaje 
para la separación entre drenes L, = 25 m 
ha permitido mostrar que el parámetro de 
escala de la radiación fractal no depende de 
la separación entre drenes. En consecuencia, 
si se tienen drenes instalados a 1.20 m 
de profundidad en un suelo de textura 
arcillosa, este parámetro puede considerarse 
proporcional a la conductividad hidráulica a 
saturación del suelo. 

Se ha mostrado que cuando las transferen- 
cias desarrolladas en la zona vadosa del 
suelo se representan de manera simplificada, 
considerando exclusivamente la ecuación 
de Boussinesq con recarga nula o igual a la 
evapotranspiración variando en el tiempo, se 
traslada la incertidumbre de lo desarrollado 
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en esta zona al valor de los parámetros que 
se estimen a partir de las pruebas de drenaje. 
Este problema puede ser minimizado con el 
modelo Richards-Boussinesq propuesto en este 
trabajo. 
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Abstract 


ZAVALA, M., SAUCEDO, H. + FUENTES, C. Fractal radiation model for agricultural 
drainage based on the Richards and Boussinesq equations. Water Technology and Sciences, 
formerly Hydraulic engineering in Mexico (in Spanish). Vol. II, No. 3, July-September, 
2011, pp. 141-157. 


The movement of water in an agricultural subsurface drainage system is analyzed using 
a model based on the coupling of the Richards equation and the Boussinesq equation, with 
variable storage capacity subject to the drain fractal boundary radiation condition. The model 
is used to describe the hydraulic operation of a drainage system installed in the field, which 
allows for studying drain separations of40 m and 25 m. Most of the hydraulic characterization 
of the system is conducted independently of drainage tests using a methodology based on 
the granulometric curve and total porosity of the soil. The experimental evolutions of the 
drained depth and the decrease in the hydraulic head associated with a 40 m separation 
between drains allow for calibrating two parameters of the proposed model: the pressure scale 
related to the soil-water retention curve and the scale parameter of the fractal radiation. The 
experimental evidence corresponding to a 25 m separation between drains enables showing 
that the scale parameter of fractal radiation is independent of the drain separation, and that 
it can be assumed to be proportional to the saturated hydraulic conductivity of the soil. It is 
shown that the classical hypotheses regarding agricultural drainage represent, in a simplified 
manner, the transfers developed in the partially-saturated zone, considering the recharge term 
in the Boussinesq equation to be null and or equal to the evapotranspiration rate; transferring 
uncertainty about what happens in the zone to the parameters estimated based on drainage 
tests. 


Keywords: vadose zone, saturated zone, shallow unconfined aquifer, recharge, discharge. 
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